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SUMMARY 

Coupled problems with various combinations of multiple physics, scales, and domains are found in 
numerous areas of science and engineering. A key challenge in the formulation and implementation 
of corresponding coupled numerical models is to facilitate the communication of information across 
physics, scale, and domain interfaces, as well as between the iterations of solvers used for response 
computations. In a probabilistic context, any information that is to be communicated between 
subproblems or iterations should be characterized by an appropriate probabilistic representation. 
Although the number of sources of uncertainty can be expected to be large in most coupled 
problems, our contention is that exchanged probabilistic information often resides in a considerably 
lower dimensional space than the sources themselves. This work thus presents an investigation into 
the characterization of the exchanged information by a reduced-dimensional representation and, in 
particular, by an adaptation of the Karhunen-Loeve decomposition. The effectiveness of the proposed 
dimension-reduction methodology is analyzed and demonstrated through a multiphysics problem 
relevant to nuclear engineering. Copyright © 2000 John Wiley & Sons, Ltd. 
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1. Introduction 

The modeling and simulation of coupled systems governed by multiple physical processes that 
may exist simultaneously across multiple scales and domains are critical tools for addressing 
numerous challenges encountered in many areas of science and engineering. However, models 
are, by definition, only approximations of their target scenarios and are thus prone to 
modeling errors. Additionally, parametric uncertainties may exist owing to various limitations 
in manufacturing and experimental methods. Uncertainty quantification (UQ) thus constitutes 
a key requirement for achieving realistic predictive simulations. 

Probability theory provides a rigorous mathematical framework for UQ, which permits 
a unified treatment of modeling errors and parametric uncertainties. The first step in a 
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probabilistic UQ analysis typically involves using methods from mathematical statistics [l|, Q 
to characterize the uncertain features associated with a model as one or more random 
variables, random fields, random matrices, or random operators. The second step is to map 
this probabilistic representation of inputs through the system model into a probabilistic 
representation of responses. This can be achieved in several ways, which include Monte Carlo 
sampling techniques Q and stochastic expansion methods. The latter typically involve the 
computation of a representation of the predictions as a polynomial chaos (PC) expansion. 
Several approaches are available to calculate the coefficients in this expansion, such as 
embedded projection 0,1^, nonintrusive projection and collocation [it-flO]. 

A key challenge in the formulation and implementation of a coupled model is to facilitate the 
communication of information across physics, scale, and domain interfaces, as well as between 
the iterations of solvers used for response computations. This information can comprise physical 
properties, energetic quantities, or solution patches, among other quantities. Although the 
number of sources of uncertainty can be expected to be large in most coupled problems, we 
believe that the exchanged information often resides in a considerably lower dimensional space 
than the sources themselves. Exchanged information can be expected to have a low effective 
stochastic dimension in multiphysics problems when this information consists of a solution 
field that has been smoothed by a forward operator, and in multiscale problems when this 
information is obtained by summarizing fine-scale quantities into a coarse-scale representation. 

In this work, we thus investigate the effectiveness of dimension-reduction techniques for 
the representation of the exchanged information. We propose to represent the exchanged 
information by an adaptation of the Karhunen-Loeve (KL) decomposition as this information 
passes from subproblem to subproblem and from iteration to iteration. When the exchanged 
information has a low effective stochastic dimension, this representation allows a reduction 
in the number of requisite stochastic degrees of freedom to be achieved while maintaining 
accuracy, thus paving the way for a solution in a reduced-dimensional s pac e, which in turn 
reduces the computational cost. It should be noted that in references Ill4l6|. the integration 
of dimension-reduction techniques in algorithms for solving stochastic partial differential 
equations has already been demonstrated; however, this paper contributes by highlighting 
the role that dimension-reduction techniques can play in solving coupled problems. 

The organization of this paper is as follows. First, in Sec. [2j we outline the proposed 
methodology. Then, in Secs.[3]and|4j we describe a version of the KL decomposition that is well- 
adapted to construct a reduced-dimensional representation of exchanged information. In Scc.0 
we focus on the implementation of our dimension-reduction methodology. Finally, in Sees. [S] 
and [3 we demonstrate the effectiveness of the proposed dimension-reduction methodology 
through an illustration problem for which we also provide numerical results. 



2. Dimension-reduction methodology 

2.1. Model problem 

This paper is devoted to the solution of a stochastic coupled model of the following form 
/(tt,a;,£) = 0, y = h{u,i), / : R r X R s ° X R m ->■ R r , h:WxR m - 



g(y,v,C) = 0, x = k(v,C), 9 ■ R r ° X ^ s x R n -> R s , k : M s x 



ro 

(1) 
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To avoid certain technicalities involved in infinite-dimensional representations, we assume that 
these equations are discretized representations of a stochastic model that couples two physics, 
two scales, two domains, or a combination of these subproblems. For instance, these equations 
may be obtained from the spatial discretization of a steady-state problem, or they may be 
the equations obtained at a single time step after the spatial and temporal discretization of 
an evolution problem. Further, we assume that the data of the first subproblem, which enter 
this subproblem as coefficients or loadings or both, depend on a finite number of uncertain 
real parameters denoted as £i,.-.,£m, and that the data of the second subproblem depend 
on a finite number of uncertain real parameters denoted as £1, . . . , £„. Lastly, we model these 
sources of uncertainty as random variables and collect them into vectors, £ = (£i,...,£ m ) 
and C = (Ci> • • • i Cn)j which are assumed to be defined on a probability triple (0, T, P) and 
considered to have values in K m and K™, respectively. A probability triple is characterized by 
its constituents: 0, a sample space of possible outcomes; T, a collection of subsets of known 
as events; and P, a probability measure. 

The stochastic coupled model ([I]) is a general bidirectionally coupled model. The solution 
variables, u, of the first subproblem, /, depend on the solution variables, v, of the second 
subproblem, g, through coupling variables, x; likewise, v depends on u through y. 

Thus, to solve this stochastic coupled model, we require to find the random variables u and v 
defined on (0,7", P) with values in W and R s such that (Q} is satisfied under the assumption 
that the stochastic coupled model is well-posed in that it admits a unique and stable solution. 

For example, if the stochastic coupled model were a fluid-structure interaction model, Eq. (JTJ) 
could represent the fluid and the structural model; the solution variables u and v could collect 
the solution fields required to describe the states of the fluid and the structure; and the coupling 
variables x and y could be the traces of the velocity field of the structure and the pressure 
field of the fluid on the deforming fluid-structure interface. 

2.2. Partitioned iterative solution 

Because a coupled model usually characterizes its response only in an implicit manner, the 
numerical solution of a coupled model typically requires an iterative method. One strategy 
could be to develop an operator-specific iterative method and associated solver for the coupled 
model as a whole. However, here, we assume that iterative methods and solvers already exist 
for each subproblem, and we therefore consider a hybrid iterative method that reuses the 
aforementioned separate solvers as steps in a global iterative method built around them to 
obtain a solution to the coupled model. The former approach is commonly labeled monolithic 
in the literature, whereas the latter is referred to as partitioned. Clearly, the key advantage of 
the partitioned method is that it enables immediate reuse of legacy software that may already 
be available to solve subproblems within their own physical and geometric domains, where one 
can account for suboperator specifics such as mesh details and problem scales. 

Let us assume that each of the aforementioned separate iterative methods is based on the 
reformulation of the associated subproblem as a fixed-point problem: 

u = a(u, x, £), y = h(u, £), a : W x R S( > x W n -> W r , h : R r x K m R r °, 
v = b(y,v,C), x = k(v,C), 6 : M r ° x R s x R™ — > M s , k : R s x R n -> R s ° . 

It should be noted that these equations can be obtained by setting a(u, v, £) = u — /(it, v, £) 
and b(u, v, £) = v — g(u, v, £), but that alternative reformulations, such as those involving a 



Copyright © 2000 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cls 



Int. J. Numer. Meth. Engng 2000; 00:0-0 



STOCHASTIC MODELING OF COUPLED PROBLEMS 



3 



direct solution of the subproblems or of their linear approximations, are often better adapted. 
We then consider the solution of the stochastic coupled model by a Gauss-Seidel iterative 
method using suitable initial values u°, v°, and x° — k(v°, £) as follows: 

v t = b( V t y-\c), **=k(v*,a 

This is not the only partitioned iterative method available; however, for simplicity, we employ 
only this method in this work. It should be noted that although we implement the proposed 
methodology using the Gauss-Seidel iterative method, one can readily use the proposed 
methodology with other iterative methods such as Jacobi, relaxation, and Newton methods. 

Further, it should be noted that the successive approximations determined by the iterative 
method ([3]) can be constructed as random variables of the following form: 

^(0)W(£(0),C(0)), 
= t/(£(0),CW); 

i.e., u e and v e can be constructed as transformations of the input random variables £ = 
• • • >£m) and C — (Ci> • • • ! Cn)' The random variables u £ and v thus exist in a solution 
space of stochastic dimension m + n. 

Finally, we acknowledge that analytical conditions under which the sequence of 
approximations generated by is ensured to converge can be established. However, a detailed 
treatment of this issue is beyond the scope of this paper, and we only note that the Banach 
contraction-mapping theorem and the attracting fixed-point theorem [17| allow one to establish 
convergence properties under global Lipschitz continuity and local differentiability conditions. 

2.3. Dimension reduction 

We believe that exchanged information often resides in a considerably lower dimensional space 
than the input sources of uncertainty themselves. Therefore, we investigate the effectiveness of 
dimension-reduction techniques for the representation of the exchanged information. Rather 
than exchanging the coupling variables x e and y and the solution variables u l and v in 
their original form, we propose to approximate these random variables by a truncated KL 
decomposition as they pass from subproblem to subproblem and from iteration to iteration. 
We specifically consider the solution of the stochastic coupled model by a Gauss-Seidel iterative 
method that accommodates a dimension-reduction technique as follows: 

u* = a(u*-^y-^,Z), y* = h{u\(-), 

where q l,d — [y e ' d ; and r e > e = [u l ~ 1 ' e - 1 x l ~ 1 ' e ] are truncated KL decompositions of 

qi = [y ;f> 1 ] and r l — [m^ , sc^ -1 ], respectively, which read as follows: 



j t (6) 

7^ = 7^ + ^^^. 



3=1 

These decompositions are described in detail in later sections. 
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It should be noted that © provides a combined reduced-dimensional representation of y l 
and v i ~ 1 in terms of a single set of reduced random variables rf — (ryf , . . . , rf^) and a combined 
reduced-dimensional representation of and in terms of a single set of reduced random 
variables t = . . . , b e ). However, this is not the only construction of a reduced-dimensional 
representation that could be considered. One can readily use the proposed methodology 
with other dimension-reduction techniques such as those involving the construction of a 
separate reduced-dimensional representation of the coupling and solution variables, with each 
representation having its own reduced random variables. 

Truncation of KL decompositions most often results in approximation errors, and we thus 
use a hat superscript to distinguish the successive approximations determined by ([5]) from 
those determined by ([3]). Further, we say that exchanged information has a low effective 
stochastic dimension when its KL decomposition ([6]) can be truncated after a few terms while 
maintaining sufficient accuracy. Random variables can be expected to have a low effective 
stochastic dimension when their components exhibit significant statistical correlation. 

Owing to the dimension reduction, the successive approximations determined by the iterative 
method ([5]) can be constructed as random variables of the following form: 

i.e., ii l can be constructed as a transformation of the input random variables £ = . . . , £ TO ) 
and the reduced random variables t = and i) e can be constructed as a 

transformation of rf = (rjf , . . . , rj l d ) and C = (Cij • • • ) Cn)- The random variable ii e thus exists 
in a space of stochastic dimension m + e and v exists in a space of stochastic dimension d + n. 

Finally, although our notations do not express a potential dependence of d and e on £, it 
should be noted that the reduced dimensions can be allowed to depend on the iteration. 

2.4. Effectiveness of the proposed dimension-reduction methodology 

The key feature of the proposed methodology is that it enables a solution of the subproblems 
in a reduced-dimensional space when the exchanged information has a low effective stochastic 
dimension. Specifically, a solution in a reduced-dimensional space is enabled when the reduced 
dimensions can be selected such that d < m and e < n while maintaining sufficient accuracy; 
refer to ^ and ([7]). This benefit is of particular significance for implementations of stochastic 
coupled models using stochastic expansion methods. These methods suffer from a curse of 
dimensionality in that their computational cost increases quickly with an increase in the 
stochastic dimension. The proposed methodology addresses the curse of dimensionality by 
mitigating the increase in stochastic dimension when information is exchanged. Although this 
prospect is the main motivation for this work, we limit ourselves in this paper to the description 
and analysis of the dimension reduction itself, and we defer the presentation of algorithms that 
exploit the dimension reduction to achieve computational gains to a later paper. 

The proposed methodology also has the potential of providing insight into the structure and 
evolution of the information that is exchanged between the subproblems and the iterations. 

Finally, it should be noted that the proposed methodology can readily be adapted to meet 
various requirements of specific applications. One could, for instance, use the KL decomposition 
to represent only those exchanged random variables that are of low effective stochastic 
dimension and adopt an alternative probabilistic representation for the remaining variables. 
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3. Karhunen-Loeve decomposition 



Here, we recall the use of the KL decomposition to construct a reduced-dimensional 
representation of a random variable q that is defined on a probability triple (0,T, P), takes 
values in a Euclidean space W" , and is of the second order: 

/ ||<?|| 2 cLP < +00, (8) 
Je 

where |j-|j denotes the Euclidean norm. 



3. 1 . Interpretation 

In the context of the proposed methodology, we think of q as a random variable that passes 
from one subproblem to another or from one iteration to the next. With reference to ©, the 
random variable q could collect the components of y and v, with w = ro + s, or the components 
of it and x, with w = r + so, at a specific iteration. 



3.2. Second-order descriptors 

The mean vector q and the covariance matrix C q of the second-order random variable q are 
defined as the w-dimensional vector and square matrix such that 



q= / qdP, (9) 
Je 

C q = f (q-q)(q-q) T dP. (10) 
3.3. Reduced- dimensional representation 

Because C q is a symmetric and positive semidefinite matrix, the solution of the eigenproblcm 

C q <f> 3 = \j<p? , (11) 

provides a set of w eigenvalues Ai > A2 > ... > X w > in addition to w 
eigenvectors dn 1 , . . . , <p w , which constitute an orthonormal basis of W" such that 

(^)V=*tf. (12) 

where is the Kronecker delta, which is equal to 1 if i = j and otherwise. The KL 
decomposition of q is then given by 

w 

Q = Q + J2\/%Vj4> J , (13) 

3=1 

where the T]j are random variables defined on (Q,T,P), with values in R, such that 

»fc = -^=(«-?)V. (14) 

V A 3 
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Further, the rjj are zero-mean and uncorrelated: 

/ VjdP = 0, 
Je 

/ Villi dP = Si 
Je 



(15) 
(16) 

The truncation of (1131) after d terms provides a reduced-dimensional representation as follows: 

d 
3=1 

in which the truncation error, because of the orthonormality properties (|12l) and (fT6)) , satisfies 

~ w 

/ \\q-q d \\ 2 dP= V (18) 

Equality (fT5)l indicates that the accuracy of (fTT)) can be improved systematically by increasing 
the number of retained terms. Further, for any orthonormal basis {e 1 , . . . ,e w } of M. w , the 
random variable q can be expanded as q = q + Y^j=i(q ~ q) T e J e-? : . The error £^ = 
S7=ci+i {q—q) T e :j e : > introduced owing to the truncation of this expansion after d terms satisfies 

T 



WXdP. 



i=d+i 



J2 (q-q) T eie>)dP 

j=d+i 



f. w 



q-qfe 3 ) dP 



j=d+i 



= ^ (^) T C q e^ 



(19) 



This expression indicates that the reduced-dimensional representation (|17p is optimal because 
from among all decompositions with d orthonormal basis vectors, this representation minimizes 
the mean-square norm of the approximation error, as shown below: 



/ \\q - q d \\ 2 dP = min / q - [ q + ~S^(q — g) T e J e- J 

ie" " {e 1 ,...,e-}6B- J& V £j 



dP. 



(20) 



(e*) T e^=(5, 



4. Proposed adaptation of the Karhunen-Loeve decomposition 

When our dimension-reduction methodology is applied to a stochastic coupled model that 
has been discretized in space and time, this methodology repeatedly requires the reduction 
of random variables with values in a finite-dimensional Euclidean space. The standard KL 
decomposition, recalled in the previous section, allows random variables with values in a finite- 
dimensional Euclidean space to be reduced in a manner that is optimal in the natural mean- 
square norm. However, it is not always suitable to use a dimension-reduction technique that 
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aims at maintaining accuracy in the natural mean-square norm; for instance, when the samples 
of the random variables to be reduced have space or time derivatives that also carry important 
pieces of information, it may be preferable to use an alternative dimension-reduction technique 
that also gives weight to maintaining accuracy in these derivatives. The significance of the 
function values, derivatives, and other hallmarks associated with the random variables to be 
reduced can be expected to be determined by the function-analytic structure that the stochastic 
model exhibited prior to its discretization. In this section, we present an adaptation of the KL 
decomposition to address this issue. Our adaptation allows a random variable that solves a 
space-time discretized stochastic model to be reduced in a manner that maintains consistency 
with the function-analytic structure exhibited by the stochastic model before discretization. 

4 ■ 1 ■ Methodology 

In the research area of stochastic modeling and analysis, the KL decomposition is known best 
as a tool for the reduction of stochastic processes that describe uncertain fields of coefficients 
and boundary conditions of stochastic partial differential equations. In such applications of the 
KL decomposition, the classical description of stochastic processes as indexed collections of 
random variables is used. In this work, we rather intend to reduce the solution of stochastic 
models or to reduce coupling variables that depend on this solution through a specific mapping. 
However, the solution of many stochastic models, including that of many stochastic partial 
differential equations, is not always amenable to a description as a classical stochastic process 
and is often better described as a random variable that takes its values in a function space. As 
we have already mentioned, the structure of this function space can be expected to determine 
the significance of the function values, the derivatives, and the various other hallmarks of the 
solution that needs to be reduced. A detailed explanation of the distinction between classical 
stochastic processes and function-space-valued random variables can be found in [18l LL9| 
and the references therein. This distinction has repercussions on the manner in which the 
second-order descriptors are defined and the properties of these second-order descriptors can 
be exploited to construct a reduced-dimensional representation. Thus, the approach adopted 
in this section is as follows. First, we describe the construction of a KL-type decomposition 
of a random variable with values in a function space. Then, we describe a construction of a 
KL-type decomposition of a random variable with values in a discrete approximation of this 
function space such that consistency with the structure of this function space is maintained. 

4-2. KL-type decomposition of function-space-valued random variables 

Let q be a random variable defined on a probability triple (0, T, P) with values in a separable 
Hilbert space H. Let (-,-)h denote the inner product on H and = \/(-,-)h denote the 
norm induced by the inner product. Let the random variable q be of the second order: 



This level of abstraction is very general: in addition to the standard Euclidean spaces of vectors 
and matrices, examples of H include spaces of square-integrable functions, spaces of sequences, 
and certain Sobolev spaces, among many other possibilities. 




(21) 
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4-2.1. Second- order descriptors. The manner in which second-order descriptors and their 
properties are defined for random variables with values in an infinite-dimensional Hilbert 
space is more complicated than the manner in which they are defined for random variables 
with values in a finite-dimensional Euclidean space. All definitions used here are consistent 
with those given in references H|-21|. Now, the mean is defined as the linear function m q 
from H into K such that 



m„(p) 



(q,p) H dP, VpeH, 



and the covariance is defined as the bilinear function c q from H x H into 



such that 



i(p,r)= ((q,p) H - m q (p))((q,r) H - m q (r))dP, Vp, r e H. 
Je 



(22) 



(23) 



We observe that these definitions are consistent with the function-analytic membership of the 
random variable q, in that the membership of the realizations of q in the Hilbert space H 
implies that q is analyzed most naturally through its projection onto fixed elements or basis 
vectors of H . Indeed, the mean function associates to any fixed element p of H the mean of the 
random variable (q,p) h obtained by projecting q onto p, and the covariance function associates 
to any fixed pair of elements p and r of H the covariance of the random variables (g,p)jj 
and (q,r)n. Clearly, the covariance function is symmetric: 

c q (p,r) =c q (r,p), Vp,reH. (24) 

Further, it is also positive: 

c q (p,p)>0, VpeH. (25) 

Equation (|2f [) and Holder's inequality indicate that m q is continuous. Thus, by Riesz's 
representation theorem, a unique vector q exists in H, i.e. the mean vector, such that 

(q,p) H = m q (p), VpeH. (26) 

Likewise, equation (I2f I) and Holder's inequality indicate that c q is continuous; and therefore, 
by Riesz's representation theorem, there exists j22| a unique linear continuous operator C q 
from H into H, i.e. the covariance operator, such that 



(C q (p),r) H = c q (p,r), Vp,reH. 



(27) 



Let {e J }^ 1 be any complete orthonormal basis of H. Because the Hilbert-Schmidt norm, 



» E ( c *( ei )> ej ) 2 H = x ( / (9-9- el )n(q-q, e j ) H dP 



(q-q, ei) 2 H dP 



dP 



\n 2 H 



(28) 



is bounded owing to (|2T|) . the covariance operator is [22[ a Hilbert-Schmidt operator. Moreover, 
because the covariance function is symmetric, the covariance operator is f8[ self-adjoint: 

(29) 



C — C l 
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where the operator C* from H into H is the adjoint operator of C q such that 

(C q (p),r) H = (p,C t q (r)) H , V Pl reH. (30) 
Because the covariance function is positive, the covariance operator is also positive: 

(C q (p),p) H > 0, VpeH. (31) 

4-2.2. Reduced- dimensional representation. Because C q is a Hilbcrt-Schmidt operator and 
thus compact and because C q is self-adjoint and positive, the solution of the eigenproblem 

C q {4?) = \j<jj (32) 
provides a denumerable set of eigenvalues m addition to a corresponding 

denumerable set of eigenmodes {(p° }SS=i; the eigenvalues are positive and square-summable 
in that X)j=i *j < +°° an( i thus Xj tends to as j goes to infinity; and the eigenmodes 
constitute a Hilbertian basis of H such that 

(<f>\^) H =6 ij . (33) 
The KL decomposition of the random variable q is then expressed as follows: 

oo 

g = g + I]VV^ ( 34 ) 

3=1 

where the r\j are random variables defined on (Q,T,P), with values in R, such that 

Vj = -j={q-Q^ j )h, (35) 
V x j 

and are zero-mean and uncorrelated, as given by (|15 p - (|16p . By the Hilbert-space orthogonal 
decomposition theorem 22j, the decomposition converges strongly: 

d 

~ dP = 0. (36) 

H 



lim / q - (q + J2 V^i 7 ?^' 



4-3. KL-type decomposition of function-space-valued random variables after discretization 

Now, let us consider a random variable q w that is still defined on the probability triple (O, T, P) 
but takes, this time, its values in a finite-dimensional subspace H w of H. Let this random 
variable q w be represented by an expansion of the following form: 

w 

q w =Y J d 1 n 3 , (37) 

3=1 

where the vectors n 1 , . . . , n w constitute a basis in H w , and let it be of the second order: 

||<f||^dP<+oo. (38) 



Let q = (gi, . . . , q w ) now be the random vector with values in M. w which collects the random 
coordinates in expansion (|37|) . In the remainder of this section, we will propose an adaptation of 
the KL decomposition which allows a reduced-dimensional representation of q to be obtained 
in a manner that is consistent with Hilbertian projections in H. 
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4-3.1. Second-order descriptors. Similarly to the KL decomposition of Sec.[3J the first step in 
the construction of our adaptation of the KL decomposition consists in determining the mean 
vector q and the covariance matrix C q of q, as defined by (|9l)- (fTUl) . 



4-3.2. Reduced- dimensional representation. Whereas the KL decomposition of Sec. |3] was 
obtained by solving the eigenproblem determined solely by the covariance matrix, our 
adaptation of the KL decomposition is rather obtained by solving the generalized eigenproblem 

W T C q Wcf> J = XjWft, (39) 

which features not only the covariance matrix C q but also the Gram matrix, 



W 



... (n\n w ) H 



(n w ,n 1 ) H ... (n w ,n w ) H 



(40) 



of the basis vectors n , . . . , n w as a weighting matrix; the Gram matrix is the w-dimensional, 
symmetric, positive definite matrix that collects the inner products of these basis vectors. 
Because C q is symmetric and positive semidefinite and W is symmetric and positive definite, 
the solution of (f3"9"| provides a set of w eigenvalues Ai > A2 > . . . > X w > in addition to w 
eigenvectors <p , . . . , <p w , which constitute, this time, a W- weighted orthonormal basis of WL W : 

WfWtf (41) 



It should be noted that the eigenvalues and eigenmodes of the generalized eigenproblem (|39|) 
usually differ from those of the eigenproblem pip , even though we use the same notations. As 
an alternative to (fTT|) . we then obtain a reduced-dimensional representation q d of q as: 

d 

q^q + Ys^Vjtf, (42) 
3=1 

where the rjj are, this time, random variables on (0,7~, P), with values in R, such that 

r h = -^=(q-q) T W<t> 3 , (43) 
V A i 

and are zero-mean and uncorrelated, as given by (H~5T) — (|T6|) . Again, it should be noted that 
the reduced random variables of the decomposition (|4"2")l usually differ from those of the 
decomposition (IT71) . even though we use the same notations. Because of the orthonormality 
properties, the truncation error incurred by q d satisfies: 



dii 2 
9-9 llw 



dP= A ^ ( 44 ) 



where y j & dP = y J e p T Wp(iP for any second-order random variable p on (Q,T,P) 

with values in M. w . It should be noted that whereas the truncation error was gauged by the 
natural mean-square norm in (|18|) . it is gauged here by the W- weighted mean-square norm. 
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Moreover, following a reasoning that is similar to the one exhibited in Sec. [3J it can readily 
be shown that q d is optimal because from among all decompositions with d orthonormal basis 
vectors, this representation minimizes the mean-square norm of the approximation error: 

2 

dP. (45) 



[ \\q - q d \\ 2 dP = mm / q - ( q + V(g - qfWe^e 



w 



(e') 1 We. 3 =Sij 

Again, it should be noted that whereas the reduced-dimensional representation (|17p 
achieved optimality in (|20[) in the natural mean-square norm, the reduced-dimensional 
representation (1421) achieves optimality here in the W-wcightcd mean-square norm. 

The interest of introducing the Gram matrix of the basis vectors is that the space of second- 
order random variables with values in W" equipped with the norm J j e \ \ ■ \ \ 2 W dP and the space 

of second-order random variables with values in H w equipped with the norm y J e 1 1 • 1 1 H dP 
share the same structure. Indeed, the function that maps any second-order random variable p 
with values in W" onto a corresponding second-order random variable p w = with 
values in H w is a structure-preserving linear bijection between these spaces such that 



^ Jjpf w dP= ] j Jjp^\\ 2 H dP, Vpe4(9,R"). (46) 

As a conclusion, when the random coordinates q = (qi, . . . , q w ) of a random variable q w = 
y]J—i qjU 3 with values in a function space H must be reduced, the use of the Gram matrix W of 
the basis vectors n 1 , . . . , n w as a weighting matrix in the construction of the KL decomposition 
allows a reduced-dimensional representation q d of q to be obtained which is consistent with 
Hilbertian projections in H and therefore optimal in a norm consistent with the norm of H . 

4-3.3. Relationship with the KL-type decomposition of function-space-valued random variables. 
Let <f" and be the mean and covariance operator of q w , defined in a manner consistent 
with the definition of second-order descriptors of function-space-valued random variables: 

(9 W ,P"'> H - / {q w , P w ) H dP, Vp w G H w , (47) 

Je 

(C™(p w ),r w ) H = [ (q w -7t,p w ) H (q w - q w ,r w ) H dP, Vp w ,r w G H w . (48) 

Owing to the relationship (|57| between the random coordinates q and the random variable q w , 
the second-order descriptors of q are related to those of q w as: 

q T Wp= {q w ,P W ) H , VpeR w , (49) 

p I W T C q Wr = (C™{p w ),r w ) H , Vp,reR™, (50) 

in which the correspondence between p and p w and between r and r w is such that p w = 
Sj=iPi n "' an d rW = S^=i r j nJ ') anc ^ ^ i s still the Gram matrix of the basis vectors. 

Further, the reduced-dimensional representation q d of the random coordinates q determines 
a corresponding reduced-dimensional representation q d > w of the random variable q w as follows: 

d 

q d > w = q w +Y,^ 1 4> J ' W , (51) 

3=1 
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where If 1 " — J2iLiQi n1 ' an d <^' ,w = Y^7=i 4>i ni - It follows from the relationship (J3HJ) between 
the mean vector of q and the mean vector of q w that q w in (f5Tj) is precisely the mean vector 
of q w . Moreover, it follows from the relationship (|5U)) between the covariance matrix of q and 
the covariance operator of q w that the Xj and the 4> 3,w in (fijTj) are precisely the eigenvalues 
and eigenmodes of the eigenproblem determined by the covariance operator of q w : 



Hence, the reduced dimensional representation q d ' w of q w deduced from the reduced- 
dimensional representation q d of q by (|5ip , is precisely the reduced-dimensional representation 
of q w that would obtained if the methodology for the construction of a KL-type decomposition 
of function-space- valued random variables of Sec 14.21 were applied to q w . 

4-4- Concluding remarks 

In this section, we described an adaptation of the KL decomposition which permits random 
variables with values in a finite-dimensional Euclidean space to be reduced in a manner that is 
optimal in a weighted mean-square norm. This adaptation is well-suited for the reduction of a 
random variable that solves a space-time discretized stochastic model, because by appropriately 
choosing the weighting matrix as the Gram matrix of the discretization basis, a reduced- 
dimensional representation is obtained which is consistent with the function-analytic structure 
that the stochastic model exhibited before its discretization. 



In this section, we provide details on the implementation of the proposed dimension-reduction 
methodology using stochastic expansion methods. 

5.1. Discretization using stochastic expansion methods 

Within the context of the model problem of Sec. HJ let {ip a , a G N" l + n } De a Hilbertian basis 
for the Hilbert space of P(££) square- integrable functions from R m +™ into R. Let this Hilbertian 
basis be indexed by multi- indices a = (ai, . . . , a m+n ) in N m+n . We will employ throughout 
this work a Hilbertian basis that is constituted of polynomials of increasing total degree, and 
we will thus refer to it as a Polynomial Chaos (PC) basis. Note that this PC basis can be 
obtained by Gram-Schmidt orthonormalization of a collection of multivariate monomials, or 
read from tables in the literature (3-fTojj if P(£.£) ^ s a "labeled" probability distribution. 

The iterative methods ([3]) and ([5]) require at each iteration that the subproblems be solved to 
obtain updated representations of the solution variables. Stochastic expansion methods involve 
the approximation of the solution variables of the subproblems by PC expansions. The PC 
basis provides the approximate representation of the solution to © as follows: 




(52) 



5. Implementation 



p 



u 




a 1=0 



(53) 



p 





a 1=0 
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and the approximate representation of the solution to (|5|) as follows: 

a|=0 

in which the summation from \a\ = to p means the summation over the finite subset of 
multi-indices a = (ai, . . . , a m+ „) in N m+ " with \a\ = ai + . . . + a m +„ < p. 

Then, the task of the solution algorithm is to compute the coordinates in these expansions, 
for which several methods are available, such as embedded projection 0, [H, nonintrusive 
projection [5|, and collocation [6 10]. We use the nonintrusive projection method in this 
work. Although we implement our dimension-reduction methodology using the nonintrusive 
projection method, it should be noted that our methodology can readily be adapted to other 
methods such as embedded projection and collocation. 

5.2. KL decomposition 

Our methodology hinges on the representation of information by a truncated KL decomposition 
as it passes from subproblem to subproblcm and from iteration to iteration. The 
implementation of this KL decomposition requires a probabilistic representation of the random 
variables to be reduced. In this work, we elected to represent random variables systematically 
by PC expansions. Here, we thus describe the implementation of the KL decomposition of a 
random variable that is represented by a PC expansion, and we show how this implementation 
naturally provides in turn a PC expansion of the reduced random variables. 

Adopting the notations used in Sees. [3] and [4] we consider the reduction of a second-order 
random variable q p that takes its values in M. w and is represented by a PC expansion of the 
form: 

<f = £ 9«V»«(€,C), (55) 

\a\=0 

We specifically consider the construction of the VF-wcighted KL decomposition presented in 
Sec. |U the standard KL decomposition presented in Sec. [3] can be recovered easily from the 
KL decomposition presented in Sec. [4] by simply setting the weighting matrix W equal to the 
identity matrix. Because of the orthonormality of the ip a , the mean q and the covariance C q 
of q p follow immediately from the PC coordinates: 

q = Qo, (56) 

C q = J2 q a ql. (57) 
|c|=l 

Then, the solution of the generalized eigenproblem W T C q W& = XjW^ provides the 
eigenvalues Ai > A2 > ■ ■ ■ \ w > and the associated eigenmodes (p 1 , . . . , <j> w required to 
construct a reduced-dimensional representation q v - d of q p as follows: 



^ = 9 + 1^. (58) 
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where the 77? are random variables with values in R such that 



v r = ^=( q r-q) T W& (59) 
V A i 



and are zero-mean and uncorrelated. By substituting (|55|) in (|59|) . a representation of each 
reduced random variable as a PC expansion is immediately obtained: 

<=E^«(tC) with r )j , a = -^=qlW<y, (60) 



|ct|=l V A J 

thus indicating that the KL decomposition of a PC expansion naturally provides a complete 
probabilistic characterization of the reduced random variables as a PC expansion. 



5. 3. Selection of the reduced dimension 

In this work, we let the solution algorithm automatically adjust the reduced dimension within 
each subproblem and at each iteration in such a way that a prescribed accuracy level is 
maintained. This iteration-dependent adjustment of the dimension reduction ensures that 
a persistent accuracy level is maintained as the successive iterations evolve. Naturally, and 
even though the truncation errors introduced owing to the dimension reductions can be 
made arbitrarily small by retaining a sufficiently large number of terms systematically, these 
truncation errors will most likely have an effect on the solution of the subproblems and will also 
propagate to the subsequent approximations generated by the iterative method. Appendix I 
provides a concise formal analysis of the effect of these truncation errors. 



6. Realization for a stochastic multiphysics problem 

We will now demonstrate the effectiveness of the proposed methodology through an illustration 
problem relevant to nuclear reactors. 

6.1. Problem formulation 

[Figure 1 about here.] 

We consider the stationary transport of neutrons in a one-dimensional reactor with 



temperature feedback |23(. Let the reactor occupy an open interval ]0, L[ (Fig.[T]). The problem 



then involves finding the temperature T and neutron flux $ such that 

"' f k^f] - h(T - r^) = -E { E f (T)<S>, 



dx \ dx . 

d*\ f \ (61) 

D(T)^j -(E a (T)- i /E f (T))$ = - S , 

under homogeneous Neumann boundary conditions. The first term on the left-hand side of the 
heat subproblem represents heat conduction, and the second term represents the transmission 
of heat to the surroundings; further, the right-hand side represents a distributed heat source 
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proportional to the neutron flux. The first term on the left-hand side of the neutronics 
subproblem represents neutron diffusion, and the second term represents the net effect of the 
absorption and generation of neutrons; further, the right-hand side represents a distributed 
neutron source. The coefficients k and h are the heat conductivity and heat transmittivity, 
respectively; the temperature is the ambient temperature; and v and Ef are the number 
of neutrons and the energy released per fission reaction, respectively. The coefficients D, E a 
and Sf are the neutron diffusion constant, fission cross section, and absorption cross section, 
respectively; these coefficients depend on the reactor temperature as follows: 



D(T(x)) = DrJ^, S a (T(x)) = S a , rc J^, E f (T(z)) = Sf.reJ^. (62) 



6.2. Deterministic weak formulation 

Let H = -ff 1 (]0, L[) be the space of functions that are sufficiently regular to describe the 
solutions. The weak formulation then involves finding T and $ in H such that 



/ k- — -dx+ h(T-T 00 )Sdx = E { X{(T)<f>Sdx, VS e H, 

J dx dx J Q J 

r L (7$ c7\T/ r L / \ r L 

/ D(T) ——dx+ (E a (T)-i/Ef(T))$*diB= / s^dx, V* e H. 
Jo dx dx Jo V / J 



(63) 



6.3. Random thermal transmittivity 

Uncertainties are incorporated by modeling the thermal transmittivity as a random 
field {h(x, ■), 1 < x < L} such that 



V j=i / 

where the £j are statistically independent uniform random variables defined on a probability 
triple (9, T, P) with values in [—1, 1] and the \/3£j are thus uniform random variables with unit 
standard deviation; further, the Xj and <fP are the eigenvalues and eigenmodes, respectively, 
of the eigenproblcm C(<^) = Xjft , where C is the covariance integral operator with 

<«) 

as the kernel; here, the parameter a is the spatial correlation length of {h(x, ■), 1 < x < L}. 
Clearly, the random field {h(x,-),l < x < L} thus obtained is such that the random 
variable h(x, •) has the mean h and coefficient of variation 6 at every position x, at least 
when the approximation error introduced owing to the truncation of the expansion after m 
terms is not taken into account. 
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6.4- Stochastic weak formulation 

The weak formulation of the stochastic problem involves finding random variables T and $ 
defined on (0,7~, P), with values in H, such that 



/ k- — -dx+ h(£){T -T ao )Sdx = E f T, f (T)$Sdx, VS e H, 
Jo dx dx J J 

f L <i$ d 1 ^ f L / \ f L 

/ D{T) ——dx+ (E a (T)-i/E f (T))$*£te= / stdx, V* e ff. 
Jo " x dx Jo \ I Jo 



(66) 



5. Discretization of space 

The finite element (FE) method is used for the discretization of space. The domain [0, L] is 
meshed using r — 1 elements of equal length. Let N\, . . . , N r then be a basis of element- wise 
linear shape functions such that Nj takes value 1 at the j-th node and at other nodes. Using 
this basis, the random temperature T and neutron flux $ are approximated as follows: 

r 

T r (x)=Y,T 1 N J (x), Tj G R, 

(67) 

* r (x) = y £$ j N J (x), &j G M. 



The FE discretization of the stochastic weak formulation (|66|) then involves finding random 
variables T = (Ti, . . . , T r ) and <fr = ($i, . . . , $ r ) defined on (6, T, P), with values in R r , which 
collect the nodal values of the random temperature and neutron flux such that 

[K + H(£)]T = q{*,T), 

[D(T) + M(T)]<f> = s. [ ' 

Here, K, H, D(T), and M(T) are r-dimensional matrices, and q($>,T) and s are r- 
dimcnsional vectors such that 

S?KS 2 = [ L k^^dx, (69) 
Jo dx dx 

SjHS 2 = f hS[S 7 2 dx, (70) 
Jo 

*jD(T)* 2 = jy(T^^dx, (71) 

*t M (T)* 2 = ^ (s a (T r ) - ^(T^l^dx, (72) 

S T q(T,<f>) = E& i {T r )<$> r S r dx + hT oa S r da i (73) 
Jo Jo 



L 



S T s= / s$ r dx, (74) 
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6.6. Reformulation as a realization of the model problem 

The aforementioned illustration problem can be reformulated as a particular realization of the 
general model problem introduced in Sec. [2] as follows: 

T = a(T,*,|), a:l r xR r xI m ^I r , 

* = b(T), b : R r -> R r , y ' 

where a(T, = [K + T) and b(T) = [D^ + M^)]"^. This reformulation 
indicates that the illustration problem is a simplified realization of the model problem, for three 
reasons. First, the data of the neutronics subproblem are not affected by their own sources 
of uncertainty £. Second, the neutronics subproblem admits a direct solution that does not 
require iteration. Lastly, the heat and neutronics subproblems are coupled directly through 
their solution variables rather than through intermediate coupling variables. 

6. 7. Discretization of the random dimension 

Because the sources of uncertainty £i , . . . , £ m are statistically independent uniform random 
variables with values in [—1,1], the PC basis {ip a , ot 6 N m } consists here of the normalized 
Legendre polynomials in m variables of increasing total degree with ipo = 1 . 

The PC basis provides approximate representations of the successive approximations 
determined by the iterative method that does not dimension reduction as follows: 

T Lv = J- zlv«(€), K G w, 

P (76) 

|c|=0 

6.8. Dimension reduction by KL decomposition of the temperature 

Now, we will demonstrate the proposed methodology by approximating the random 
temperature by a truncated KL decomposition as it is communicated from the heat to the 
neutronics subproblem. 

The PC basis provides approximate representations of the successive approximations 
determined by the iterative method involving dimension reduction as follows: 

|c|=0 

The mean and covariance of T e - p are given by: 

T e = T e , (78) 

p 



E TUTiV- (79) 
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Further, let the r-dimensional square matrix W be the Gram matrix of the FE basis, i.e., 

'(N^N^h ... (N u N r ) H -] 



W = 



(N r ,N r ) H 



(80) 



where the inner product (•, -)h is such that (Si, S2)h = Jq SxS^dx + (dSi / dx)(dS2/ dx)dx 
for any pair S\ and S% of functions in H . The solution of the generalized eigenproblem 
W T C e T W<j) : '- e — \ l ^W^P' 1 then provides the eigenvalues A^ and the associated eigenmodes cfi^ 

required to construct a reduced-dimensional representation T l,v,d of T l ' p as follows: 



n£,p,d 



d 
3=1 



where the rf^ p are random variables defined on (0, T, P), with values in R, such that 



e,p 



1 'T^-T^) T W^' 



By substituting (|77p in (1521 . a representation of the ?7^' p as a PC expansion is obtained: 



V"' p = E ^^a«) with rf c 



(81) 



(82) 



(83) 



thus completely characterizing the reduced random variables as a PC expansion. 

It should be noted that the random neutron flux, in principle, could also be reduced as it 
passes from the neutronics subproblem to the heat subproblem. However, because the data of 
the neutronics subproblem are not affected by their own sources of uncertainty, a reduction of 
the random neutron flux would not lower the number of sources of uncertainty that enter the 
heat subproblem, and thus would not pave the way for a solution of the heat subproblem in a 
reduced-dimensional space. This extension is therefore not demonstrated. 

6.9. Selection of the reduced dimension 

At each iteration, we select the number of terms retained in (f5Tj) by the KL 
decomposition T l ' p ' d of T l ' p as the smallest dimension d that satisfies the following condition: 



dP < toL 



w 



WeN, 



(84) 



where tol is a prescribed tolerance level. Clearly, this criterion may result in the dependence 
of the reduced dimension d on the iteration I. 



6.10. Concluding remarks 

Algorithm[T]outlines an implementation wherein the nonintrusive projection method is adopted 
for the discretization of the random dimension. This algorithm requires a quadrature rule for 
integration with respect to the multidimensional uniform probability distribution of the input 
random variables. In this work, we use a sparse-grid Gauss-Legendre quadrature rule 24fl . 
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Input : KL truncation tolerance tol; 

PC basis {ipa, < \a\ < p} up to total degree p w.r.t. P^; 
Quadrature rule {(£*., w k ), 1 < k < u} of level p + 1 w.r.t. P^; 

1 = 1: 
repeat 

heat subproblem 
for k = 1 to v do 



Solve 



K 



end 



+ ff(€ fc )]^(« fc ) = 9(2 w - 1 " , « fc ).*'- 1 " , « fc )); 

Compute PC coordinates of T e ' p by nonintrusive projection: 
end 

dimension reduction 

Compute mean T e = f e and covariance C e ~ = Ef a |=i T , i(T'i) T ; 
Solve cigenproblcm W^C^W^ 3 ' 1 = X^Wtf'*; 
Choose d such that £f Q , =1 (Ti) T WTl - E^=i A j < 
Compute PC coordinates of r^' p by = (T^^W^ 5 '* for j = 1 to d; 
end 

neutronics subproblem 
for k = 1 to do 

Solve 



£>(T^' d (£ fc )) + M(T^ d (U)\* e (t; k ) = «, 

with t^"(U = T e + T:U >/^(Efa|=i^«« fc ))^; 

end 

Compute PC coordinates of <& £,p by nonintrusive projection: 



end 

£ = £+1; 
until (convergence); 



X)**«*)iM€fc)«'*; 
fc=i 



Algorithm 1: Implementation of the illustration problem. 



7. Numerical results 

We obtained results using the following properties. We assumed the reactor to have a 
length of L = 100 [cm]. Further, we assumed a deterministic and position-independent 
neutron-diffusion constant D re f = 2.2 [cm]; absorption cross section E a>re f = 0.0195 [cm -1 ]; 
fission cross section £ a ,ref = 0.0075 [cm -1 ]; multiplication factor v = 2.2; neutron 
source s = 5.0-E11 [neutrons/s/cm ]; ambient temperature = 390 [K]; fission energy Ef = 
3.0E-U [J/neutrons]; and temperatures T rof = 390 [K], T min = 390 [K], and T max = 1000 [K]. 
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[Figure 2 about here.] 

In addition, we used a thermal transmittivity random field with position-independent 
mean h = 0.17 [J/K/cm 3 /s], spatial correlation length a = 15 [cm], and coefficient of 
variation 5 = 10%. We retained m = 10 terms in expansion (|64l) . Figure (5Ja) shows a few 
sample paths of the random field {h(x, •), < x < L} thus obtained. Figure HJ^b) shows the 
10 largest magnitude eigenvalues of the covariance integral operator. 

7.1. Monte Carlo sampling implementation 

[Figure 3 about here.] 

First, we carried out a Monte Carlo simulation. We generated MC = 100, 000 sample 
paths of the thermal transmittivity random field. Then, for each of these sample paths, we 
constructed the associated deterministic multiphysics model, each of which we solved using the 
FE method for the spatial discretization and Gauss-Seidel iteration as the iterative method. 
We systematically obtained converged results for r — 1 = 40 finite elements and 20 iterations. 

Figure [3] shows a few samples of the random temperature and neutron flux thus obtained for 
the values of k = 100 [J/K/cm/s] and k = 1 [J/K/cm/s]. This figure shows that the samples 
of the random temperature became less smooth as the value of the thermal conductivity was 
decreased; i.e., the samples exhibited more rapid oscillations with respect to the position in 
the reactor as the significance of the diffusion of heat was decreased. 

7.2. PC-based implementation not involving dimension reduction 

Next, we implemented a PC-based iterative method that did not involve dimension reduction. 
We adopted the nonintrusive stochastic projection method for the discretization of the random 
dimension. It should be noted that this implementation amounts to the implementation 
described in Algorithm [TJ provided that the dimension-reduction step is not carried out, or, 
equivalently, that the reduced dimension d is chosen equal to r at each iteration. 

[Figure 4 about here.] 

[Figure 5 about here.] 

[Figure 6 about here.] 

Figure 0] shows the convergence of the iterative method as a function of the number of 
iterations. The iterative method converged at a linear rate up to approximately iteration I = 10 
for k = 100 [J/K/cm/s] and iteration I = 15 for k = 1 [J/K/cm/s], after which linear-solver 
tolerances became dominant and prevented further convergence. 

Figures [5] and E] show the convergence of the solution as a function of the total degree at which 
the PC expansions are truncated; note that the superscript oo is used in the figure captions 
to indicate convergence with respect to the number of iterations. The distance between the 
solutions obtained through the Monte Carlo and the PC-based simulation that did not involve 
dimension reduction decreased at an exponential rate as the total degree was increased. 
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7.3. PC-based implementation involving dimension reduction 

Lastly, we implemented the PC-based iterative method involving dimension reduction. This 
implementation corresponded exactly to Algorithm [T] We obtained the results to follow using 
PC expansions truncated at total degree p = 4 and, with reference to ([541 . using a range 
of values for the error tolerance level tol adopted to determine the reduced dimension at 
each iteration. We discuss convergence as a function of the error tolerance level later. Now, 



we present detailed results obtained for tol — 0.90 x a\, where or = yJ2\ a \=i ll^allw 

is the mean-square norm of the fluctuating part of the random temperature obtained by 
the PC-based iterative method that did not involve dimension reduction: specifically, o~t = 



132.54 [JK/cm 2 /s] for k = 100 [J/K/cm/s] and a T = 201.18 [JK/cm 2 /s] for fc = 1 [J/K/cm/s]. 



Figure [7] shows the convergence of the iterative method as a function of the number 
of iterations. As in the case shown in Fig. 21 the iterative method converged at a linear 
rate up to approximately iteration £ = 10 for k = 100 [J/K/cm/s] and iteration I = 15 
for k = 1 [J/K/cm/s]. 



Figures [8] and [9] show a few components of the KL decomposition of the random temperature 
obtained at the last iteration, namely, the mean temperature, the 10 largest magnitude 
eigenvalues, and the eigenmodes and reduced random variables associated with the two largest 
magnitude eigenvalues. We can observe that the eigenvalues obtained for k = 100 [J/K/cm/s] 
(Fig. [He)) decay at a faster rate than those obtained for k = 1 [J/K/cm/s] (Fig. HDJd)). This 
behavior of the eigenvalue decay rate is consistent with our earlier observation that the samples 
of the random temperature become less smooth and therefore exhibits less spatial correlation 
as the thermal conductivity is decreased. 



Figure [10] shows a few samples of the random temperature and neutron flux deduced from 
the PC expansions obtained as the output of the solution algorithm. The samples of the input 
random variables used to synthesize the samples of the random temperature and neutron flux 
shown in Fig. [TU] were identical to those used to generate the samples shown in Fig. [31 The 
similarity of the samples in Figs. [3] and [TT] indicates that the PC-based surrogate model not 
only provides an accurate global representation of the multiphysics model but is also capable 
of accurately reproducing a sample-wise response. 

7.4- Convergence analysis 




[Figure 7 about here.] 



[Figure 8 about here.] 



[Figure 9 about here.] 



[Figure 10 about here.] 



[Figure 11 about here.] 



[Figure 12 about here. 
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We repeated the PC-based simulation involving dimension reduction for the error tolerance 
levels tol = 0.90 x erf., tol — 0.95 x erf*, and tol — 0.99 x erf,; each of these levels corresponded to 
an increased accuracy that the KL decomposition was required to maintain at each iteration. 
Figure [TT] indicates that the KL decomposition systematically retained more terms when 
the error tolerance level was set to a lower value, and, thus, higher accuracy was required. 
Moreover, the comparison of Figs. ITTT a) and Efb) reveals that the reduced-dimensional 
representations of the random temperature retained less terms for k = 100 [J/K/cm/s] than 
for k — 1 [J/K/cm/s], which is consistent with our earlier observation that the eigenvalues 
of the KL decomposition obtained for k = 100 [J/K/cm/s] decay at a faster rate than those 
obtained for k = 1 [J/K/cm/s]. 

Figure [12] shows that the distance between the successive approximations obtained by the 
iterative method that involved dimension reduction and the one that did not involve dimension 
reduction remained bounded as the iterations progressed, and that this distance can be reduced 
by improving the accuracy of the KL decomposition by decreasing the error tolerance level. 



7.5. Concluding remarks 

In this section, we demonstrated the effectiveness of the proposed dimension-reduction 
methodology through a stochastic multiphysics model of the transport of neutrons in a reactor 
with temperature feedback. A study of the parameters of the model indicated that the random 
temperature exchanged between the neutronics and the heat subproblem was of a low effective 
stochastic dimension when it was smoothed by the effect of significant heat diffusion terms. 
Numerical results showed that the KL decomposition could then extract a low-dimensional 
representation of the random temperature as it passed between the subproblems while accuracy 
was maintained. The convergence study showed that the accuracy of the solution could be 
increased systematically by retaining more terms in the KL decomposition. 



8. Conclusion 

While most coupled models can be expected to be affected by a large number of sources 
of uncertainty, information exchanged between subproblems and successive iterations often 
resides in a considerably lower dimensional space than the sources themselves. This argument 
calls for the use of dimension-reduction techniques for the representation of exchanged 
information. In this work, we described the formulation and implementation of an adaptation 
of the KL decomposition to represent information as it passes from subproblem to subproblem 
and from iteration to iteration during the numerical simulation of a stochastic coupled model. 

Implementations can capitalize on the reduced-dimensional interface created between 
the subproblems and iterations to enable a more computationally efficient solution of the 
subproblems in a reduced-dimensional space. This can be achieved by performing key 
algorithmic operations in terms of the reduced stochastic degrees of freedom of the exchanged 
uncertainty representations. This is the main computational benefit of the proposed dimension- 
reduction methodology, and we plan to explore solution algorithms of this form in future work. 
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APPENDIX 



The approximation of a random variable by a truncated KL decomposition as it is exchanged 
most often results in a truncation error, which will most likely have an effect on the solution of 
the subproblems and will also propagate to the subsequent approximations generated by the 
iterative method. This appendix concisely analyzes the effect of this truncation error. 

1.1. Objective 

The probability theory offers several ways in which two random variables can be considered 
to be equal or a sequence of random variables can be considered to converge, namely the 
almost sure, mean-square, and distributional senses, among other ways. An error analysis can 
be conducted usefully from any of these perspectives. We adopt a mean-square setting because 
it perfectly agrees with the convergence properties of the KL decomposition. Specifically, 
in a mean-square sense, we examine the distance between the sequences {{u l , v >)}'^, 1 
and {(u e , i) determined by the iterative method that does not involve dimension direction 

and by the one that does involve dimension reduction, respectively. 

1.2. Contraction-type assumption 

We refer to the mapping that transforms any random variable (it, v,y,x) defined on (0, 7", P) 
with values in l r xR s xK r ° xl s " into a random variable (a(u,x, £),b(y, v,£),h(u,g),k(v,Q) 
defined on (0, T, P) with values in l r xl s x M. r ° x M s ° as the iteration mapping associated with 
the iterative method ([3|) that does not involve dimension reduction. Clearly, an error analysis 
requires assumptions concerning the sensitivity of this mapping with respect to its arguments. 
In this analysis, we assume the iteration mapping to be a contraction mapping. Specifically, 
we assume that the iteration mapping is a contraction mapping in the block-maximum norm, 
in that there exists a Lipschitz continuity modulus a in [0, 1) such that 



max(||a(u, x, £) -a(u, x , £)\\ r ,\\b{y, v, Q -b(y', v, Q\\ s ,\\h(u, £) -h(u, £) || ,||fc(«, C) -*(«', C)|L) 



for all the second-order random variables (u,v,y,x) and (u' , v', y', x') defined on (0,7", P) 
with values in W x R s x R r ° x R s °; here, |j-|| r , ||-|| s , ||-|| , and ||-|| are suitably weighted 
mean-square norms on the spaces of second-order random variables on (0, T, P) valued in R r , 
K s , R r °, and R s °. This assumption is justified in the current context by the fact that this 
assumption is consistent with conditions that ensure the mean-square convergence of the 
sequence generated by the iterative method ([3]) that does not involve dimension reduction: by 
the Banach contraction mapping theorem, a sufficient condition for {(u e ,v^'}^_ 1 to converge 
in mean square is that the iteration mapping should map a non-empty closed set of random 
variables into itself and that (|85|l should be fulfilled. 

1.3. Adaptive selection of the reduced dimension and induced error bound 

If the iterative method ((3|) not involving dimension reduction and the iterative method ([5]) 
involving dimension reduction are initialized identically, if the iteration mapping is a 
contraction mapping in the block-maximum norm such that (|85p is fulfilled with continuity 
modulus a in [0, 1), and if the reduced dimensions, denoted here as dg and ei, are selected at 
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each iteration so as to achieve systematically a prescribed accuracy tol > such that 
maxfllw^ 1 -ft*- 1 ' 8 *!! J^- 1 -^- 1 '^!! )<tol, WeN, 

Ml Mr ' II 1 1 s ' — ' ' /ao\ 

max(\\y e -y e ' d '\\ ro , llv 1 ' 1 -v*- 1 ' 4 *^) < tol, WeN, 

then the distance between the sequences of successive approximations generated by the iterative 
method that does not involve dimension reduction and the one that does involve dimension 
reduction satisfies 

max(||i/-i/|| , M-v'll ) < -^-tol, WeN. (87) 

" Mr II 1 — a 

It should be emphasized that the reduced dimensions dg and eg can always be chosen such 
that (|86|) is fulfilled owing to the key properties (fT8f and (|44|) of the KL decomposition. 



1.4- Proof of the error bound 

From the definition of the iteration processes © and ([5]), from the contraction- type 
assumption ([551) . and from the triangle inequality, it is observed that 

max(||t/ - u\, \\v e - 6*11 J < amaxfllti' -1 - u' _1 || r> - * £_1 || s ) + 2atol. (88) 

If (u°,v°) = (u°,v°), the repeated application of this inequality yields 

max(||i/-« £ || r ,||t/-^|| s ) <J22a e ~ k+1 tol, (89) 

fc=i 

which subsequently yields (|87l) because by the convergence of geometric series for < a < 1 , 
we obtain \im.e^ +00 $^t_i a l ~ k+1 — a/ (1 — a). This reasoning indicates that a contraction-type 
assumption is sufficient to obtain a result of the form given by (|87j) because it ensures that the 
effect of any approximation error introduced owing to the truncation of a KL decomposition is 
systematically diminished, and not amplified, as the error propagates to subsequent iterations. 

1.5. Concluding remarks 

This appendix showed that under a contraction-type assumption, the sequence of successive 
approximations determined by the iterative method involving dimension reduction can be 
brought as close as desired to the sequence determined by the iterative method not involving 
dimension reduction by suitably adjusting the reduced dimension at each iteration. 

Although this appendix provides a useful result, it is not a comprehensive convergence 
and error analysis of the proposed dimension-reduction methodology. Further research can be 
conducted to analyze convergence properties and error estimates under weakened assumptions 
such as those involving nonuniform contraction-type or local differentiability properties, for 



instance, by following the approaches given in 25-27]. Moreover, a posteriori error estimators 



can be developed, for instance, by following the approach given in [25 
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Figure 1. Schematic representation of the problem. 
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(b) Eigenvalues. 

Figure 2. Thermal transmittivity random field: (a) five samples and (b) ten largest magnitude 

eigenvalues of the covariance integral operator. 
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(c) Neutron flux for k = 100 [J/K/cm/s]. (d) Neutron flux for k = 1 [J/K/cm/s]. 

Figure 3. Monte Carlo simulation: five samples of the solution. 
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Number of iterations [-] Number of iterations [-] 



(a) Convergence for k = 100 [J/K/cm/s]. (b) Convergence for k = 1 [J/K/cm/s]. 
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* -> \/Ef aN o - ^fw/^U^o (squares). 

Figure 4. Simulation not involving dimension reduction: convergence with respect to the number of 

iterations. 
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(a) Convergence for k = 100 [J/K/cm/s]. (b) Convergence for k = 1 [J/K/cm/s]. 
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Figure 5. Simulation not involving dimension reduction: convergence with respect to the polynomial 

order. 
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Figure 6. Simulation not involving dimension reduction: convergence of the error estimates with respect 

to the number of samples for p = 4. 
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Figure 7. Simulation involving dimension reduction: convergence with respect to the number of 

iterations. 
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(c) Eigenvalues for k = 100 [J/K/cm/s]. 
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(d) Eigenvalues for k = 1 [J/K/cm/s]. 
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Figure 8. Simulation involving dimension reduction: mean, ten largest magnitude eigenvalues, and 
first eigenmode of the KL decomposition of the random temperature. 
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(a) Second eigenmode for k = 100 [J/K/cm/s], 
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(b) Second eigenmode for k = 1 [J/K/cm/s]. 
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(c) First reduced r.v. for fe = 100 [J/K/cm/s]. 
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(d) First reduced r.v. for k = 1 [J/K/cm/s]. 
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Figure 9. Simulation involving dimension reduction: second eigenmode and PC coordinates of the first 
and second reduced random variables of the KL decomposition of the random temperature. 
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(c) Neutron flux for k = 100 [J/K/cm/s]. (d) Neutron flux for k = 1 [J/K/cm/s]. 

Figure 10. Simulation involving dimension reduction: five samples of the solution. 



Copyright © 2000 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cls 



Int. J. Numer. Meth. Engng 2000; 00:0-0 



38 



FIGURES 



•a 6- 
I 

■8 4 



♦ ♦ ♦ ♦ ♦♦♦♦♦♦♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ 



♦ ♦♦♦♦♦♦♦♦♦ 
♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ 



a 



o 



10 



15 



Number of iterations [-J 

(a) Reduced dimension for k = 100 [J/K/cm/s]. 



20 











10 



15 



Number of iterations [-J 

(b) Reduced dimension for k = 1 [J/K/cm/s]. 



20 



Figure 11. Convergence analysis: reduced dimension as a function of the iteration for tol = 0.90 x a\ 
(circles), tol — 0.95 x a% (squares), and tol = 0.99 x o% (diamonds). 
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Figure 12. Convergence analysis: mean-square distance between the successive approximations 
generated by the simulations involving and not involving dimension reduction for tol — 0.90 x a? 
(circles), tol — 0.95 x o\ (squares), and tol = 0.99 x (diamonds). 
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